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1 Introduction 



Many first order phase transitions in Nature happen so that the external parameter 
which drives the transition varies on timescales several orders of magnitude longer 
than the microscopic interaction timescale. Examples include freezing or vaporization 
of liquids, several magnetization transitions, and a number of phase transitions in 
the early Universe. Despite the very slow change in the control parameter, these 
transitions do not happen immediately when the thermodynamic transition point of 
the control parameter is reached, but the system can remain in the old, metastable 
phase for an extremely long time compared to the interaction timescales. (For the case 
where a temperature change induces the phase transition, this is called supercooling or 
superheating.) In homogeneous systems, the transition happens through spontaneous 
nucleation of droplets (or bubbles) larger than a critical size, and subsequent droplet 
growth and coalescence. 

Let us, for a moment, specialize to temperature driven phase transitions. Given a 
fixed cooling rate, the depth of supercooling is determined by the temperature depen- 
dent nucleation rate of critical droplets. The more supercooling, the farther the system 
is driven away from thermodynamical equilibrium, and when the transition finally hap- 
pens, the more dramatic the consequences will be: hydrodynamic flows, shock waves, 
entropy production, generation of new length scales etc. In order to make quantitative 
predictions we need to know the nucleation rate accurately.^ 

A statistical theory of droplet nucleation in liquids was formulated already in 1935 by 
Becker and Doring 0]. Later, Cahn and Hilliard gave a phenomenological descrip- 
tion of the nucleation in Ginzburg-Landau theory. Finally, the theory of nucleation 
was given a firm field theory foundation by Langer [^] (for a review, see Ref. 0). 
Langer's method was generalized to relativistic quantum field theories by Callan and 
Coleman P and Voloshin [0 , and extended to finite temperature quantum field theories 
by Affleck § and Linde [|. 

A fully consistent application of Langer's theory can be difflcult. This is especially so 
in theories with a radiatively induced first order transition. The tree-level Lagrangian 
of these theories does not show a first order transition; the transition becomes apparent 
— the effective potential develops more than one minimum — only after some field 
fluctuations have been integrated over. On the other hand, in nucleation theory we start 
from a potential which has two minima, one of them metastable, and flnd the classical 

* As an example, the electroweak phase transition in the early Universe may have been of first order, 
in which case the dominance of matter over antimatter may arise as a result of strong metastability. 
However, the efficacy of the generation mechanism depends sensitively on the quantitative details of 
the transition and the degree of the supercooling. For a review, see, for example . 
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droplet solution of the equations of motion, with energy H . The nucleation rate (per 
time and volume) is now proportional to the Boltzmann factor: F = D exp{—H/T). 
Langer's theory gives us a recipe how to evaluate the constant of proportionality D; an 
essential ingredient is the determinant of the fluctuations around the classical critical 
droplet solution. The problem with the radiatively induced transitions is now apparent: 
we already "used up" some fluctuations for evaluating the effective potential. Thus, 
care has to be taken to ensure the proper counting of the fluctuations. Note also that 
the fluctuation determinant represents only a one loop treatment of fluctuations and 
ignores nonlinearities between fluctuations; in some cases these are very important, 
and a high loop or nonperturbative treatment is needed. 

In this paper we study the nucleation rate in the cubic anisotropy model, a theory 
with two scalar fields, using numerical lattice simulations. It is one of the simplest field 
theories with a radiatively generated first order phase transition, and thus well suited 
for both analytical and numerical studies. Partial results have already been published 
in Ref. |T^. The phase structure of the theory has been studied using e expansion 
techniques [0, |12| and Monte Carlo simulations |jl3|. The droplet nucleation rate itself 
has been studied using coarse grained effective potentials . 



In principle, it is relatively straightforward to study the nucleation rate using stan- 
dard real-time lattice simulations, provided that the metastability is not strong: one 
simply prepares initial configurations in the metastable phase, lets them evolve in time 
and waits for the transitions to happen. This method has been used in simulations 



of kinetic discrete spin models |16l and in scalar field theories • However, it is 



applicable only if the metastability time scale is, at most, a few orders of magnitude 
longer than the microscopic interaction time scale. Otherwise the system simply re- 
mains metastable for the duration of any practical computer simulation. This is the 
situation for many phase transitions in Nature. 

In contrast, the method described in this paper is applicable to first order transitions 
with almost arbitrarily strong metastability.^ Our method is by no means specific to 
the cubic anisotropy model; it can be used to study the nucleation rate in any theory, 
provided that: 

(1) The thermodynamics of the system is well described by an (essentially) classical 
partition function, and both the thermodynamics and the real time evolution of the 
system can be reliably simulated on the lattice; 

(2) there is an "order parameter like" observable which can resolve the phases and 

^ We note that algorithms exist in the hterature for studying metastable states in some discrete 
spin models (for example, the Ising model), and for specific choices for update dynamics. For a review, 
see Rcf. ITHl and references therein. 
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the potential barrier between them sufficiently accurately; and 
(3) the barrier is large and the nucleation rate is small. 

In a simultaneous calculation this method was applied to SU(2) gauge + Higgs 
theory, which is an effective theory for the electroweak phase transition [0. With 
suitable generalizations one can apply the method also for almost any strongly sup- 
pressed process in various theories on and off the lattice. Indeed, the method here is 
closely related to the one used for studying the "broken phase sphaleron rate" in the 
electroweak theory ||2CI|| . 

The Monte Carlo simulations are free of the ambiguities and problems which analyti- 
cal calculations suffer from in practice. Relatively modest computational resources give 
us the nucleation rate, as a function of the supercooling, quite accurately — the cal- 
culations in this work were performed using only standard workstations and PC's. We 
compare the results with various (semi) analytical approaches: the thin wall approxi- 
mation, where we use latent heat and surface tension determined by lattice simulations; 
purely perturbative analysis; and the coarse grained effective potential calculation by 
Strumia and Tetradis |T^. The results from these calculations show huge sensitivity 
to various approximations in analytical calculations. On the other hand, we do not 
observe any signs of the breakdown of hanger's theory of homogeneous nucleation in 
radiatively generated transitions, reported in Ref. ||14|| . 

This paper is structured as follows: in Sec. ^ we give a detailed description of our 
lattice approach. Since the method is quite general, the discussion is completely in- 
dependent of the specifics of the model studied and can be read independently of the 
rest of the paper. This section is a reformulation and development of the discussion in 
Sec. II in Ref. The cubic anisotropy model, its lattice discretization, and its real 
time evolution is presented in Sec. |^. The thermodynamics of the model are studied 
nonperturbatively on the lattice in Sec. ^, and the nucleation rate is determined in 
Sec. ^. The comparison with analytical and semianalytical approaches is in Sec. ^, and 
finally, we conclude in Sec. 0. 



2 How to calculate the droplet nucleation rate 

In this section we describe the general framework of the method for calculating the 
droplet nucleation rate with lattice Monte Carlo methods. The calculation can be split 
into two parts: (1) the measurement of the probability of the critical droplets, and (2) 
the calculation of the dynamical factors which convert the probability into a rate. In 
the case we are interested in, very weak supercooling, the probabilistic suppression is 
by far the dominant part. Thus, for many purposes, a sufficiently accurate answer can 
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be obtained by calculating the droplet probability in the way described below, and 
multiplying it with suitably chosen dimensional factors in order to convert the proba- 
bility into a rate. However, without the correctly calculated dynamical information the 
resulting rate can be off by several orders of magnitude. Moreover, only by including 
the dynamical information does the result become truly independent of the choice of 
order parameter and other details of the procedure. 

In order to make the discussion more concrete, we shall adopt the terminology from 
the cubic anisotropy model (described in Sec. |^), and consider a first order symmetry 
breaking transition from a symmetric phase to a broken phase, which occurs when 
a control parameter is lowered below the (adiabatic, infinite volume) transition 
value m^. However, we emphasize that the method described here is not specific to 
the cubic anisotropy model or even to order-disorder transitions in general; indeed, by 
suitable generalization it can be applied to any system where we have two or more 
states separated by a large potential barrier, and where the system evolves in phase 
space with some dynamical prescription amenable to computer simulations. 



2.1 The probability of a critical droplet 

Let us consider a system in a finite cubical volume V = L^, with periodic boundary 
conditions. At the transition point = the probability distribution P{6) of a 
suitably chosen volume averaged "order parameter"^ 

= ^1 d^xe{x), (2.1) 

attains the familiar two-peak form, see Fig. 0. The two peaks correspond to the sjth- 
metric and broken phases, with order parameter expectation values 6'symm. and ^broken, 
respectively. The intermediate region between the peaks consists of configurations 
with co-existing domains of the symmetric and the broken phases. We assume here 
that the system size is much larger than the thickness of the phase interfaces. The 
probability distribution can also be written in terms of the constrained free energy : 
P{0) oc exp[— F(^^)]. (Here and throughout we rescale free energies F and energies H 
by T so they are dimensionless.) 

^What we mean here by the "order parameter" is a quantity which has clearly different expectation 
values in the two phases, but by no means do we require that it is a true mathematical order parameter 
in the sense that it would have a zero expectation value in one phase and a non-zero value in the 
other. Thus, the transition here can be one without global symmetry breaking (for example, a liquid- 
gas transition where the density can act as an order parameter). The choice of the order parameter 



is discussed in more detail in Sec. 4.1 
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The shape of P{0) can be readily understood by considering the following "thin- 
wall" mean-field approximation: the order parameter in the bulk symmetric and bro- 
ken phases is fixed to respective expectation values, and the phase interfaces have 
zero thickness, with a surface tension a. The free energy densities of the symmetric 
and broken phases are equal, whereas the mixed phase has an additional free energy 
contribution, which equals a times the area of the phase interfaces. The geometry of 
the mixed symmetric and broken phase domains now settles in configurations which 
minimize the surface area. In this case the value of the normalized order parameter 
X = {9 — 6'symm.)/(6'broken " ^symm.) cxactly cquals the volume fraction filled by the 
broken phase. 

Given a value of X, it is a straightforward exercise to find the mixed phase geometry 
which minimizes the area. When X is small, we have a small, spherical droplet of 
the broken phase embedded in the bulk of the symmetric phase. However, when the 
diameter of the droplet approaches the length of the box L, it becomes advantageous 
for the broken phase to form a cylinder spanning the box along one direction — because 
of the periodic boundary conditions, the cylinder does not have any "end caps," and its 
area is smaller for a given volume. Analogously, when X is further increased, the broken 
phase forms a slab which spans the volume along two directions. When X > 1/2, 
the same sequence occurs again, with the symmetric and broken phases interchanged. 
Quantitatively, the stability ranges for different mixed phase geometries, in the thin 
wall approximation, follows: 

(1) droplet: < X < f , area = {SQnX^y/^L^ 

(2) cylinder: |f < X < i, area = 2{7rXy/^L^ 

(3) slab: ^ < X < 1 - i, area = 2^^ . 

These are shown with the thin lines in Fig. |I|. Note that the change in geometry is 
by no means continuous: for instance, the change from a droplet into a cylinder at 
X = 47r/81 must occur via a deformation which increases the area substantially. This 
causes an energy barrier between different mixed phase topologies, and the topology 
sectors have overlapping metastability branches.]] 

What happens to the probability distribution when we lower slightly from the 
transition value m^? Now the free energy densities in the symmetric and broken phases 

^Indeed, multicanonical Monte Carlo simulations (to be discussed in Sec. |4.2| ), which, by con- 
struction, overcome the huge probability suppression of the mixed phase, still suffer from the milder 
tunneling barriers associated with the changes of topology. This is because the order parameter 
remains constant when the topology change occurs. Thus, in multicanonical simulations at large vol- 
umes there is a clear reluctance for the system to go from the droplet branch into the slab branch, for 
example. 
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Figure 1: Order parameter probability distribution P{6) at the transition point 
in a cubic box with periodic boundary conditions. Thick hne: A schematic plot of 
the distribution on a large but finite volume. Thin line: "thin wall" approximation, 
consisting of droplet, cylinder and slab topologies. 

are different: 

= 6mH, (2.2) 

where Sm'^ = — and i is proportional to the "latent heat" of the transition 
(remember that is the temperature parameter). Remaining in our simple thin wall 
approximation, and further assuming that a remains constant as we change m^, the free 
energy of the mixed phase configurations has an additional contribution proportional 
to the volume: 

F(X, m^) = F{X, ml) + 6m^iXL^ + const. (2.3) 

In particular, the free energy of a broken phase spherical droplet of volume Vd and area 
Ad becomes 

Fd™piet(X, m^) - F„, = ~5mHVd + a Ad = -6m^iXL' + a(367r) 1/2^2/3^2 (2.4) 

The droplet free energy has a maximum Fc = IQna^ / {3{6m'^)'^£'^) at droplet radius 
Rd = 2a / {—5m?€). This is the critical droplet size and free energy at supercooling 5m^: 
droplets smaller than this size can continuously reduce their free energy by shrinking, 
and larger droplets can by growing. The droplet free energy is shown schematically in 
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Figure 2: The constrained free energy F{X) 
5m?. 



logP(X) at a small supercooling 



Fig. 1^. In a finite box the droplet volume should not exceed the droplet ^ cylinder 
stability limit L^47r/81 ~ 0.155L^, as otherwise the critical "droplet" is a cylinder 
and the box size is essential to determining its behavior, which is unphysical if we are 
interested in the behavior in the thermodynamic limit. Applied to the critical droplet, 
this implies that the inequality 5m?'^ > Qa / L should be satisfied. 

Naturally, in a more realistic case the situation is not as simple as in the thin wall 
case. The phase interfaces have a finite thickness, the shape of the mixed phase domains 
fluctuates, as also does the value of the order parameter, averaged over the volumes of 
the pure phase domains. Nevertheless, provided that the volume is sufficiently large, 
the basic topological structure outlined above remains — even if the shape of the 
cylinders, droplets and slabs becomes distorted, the magnitude of the fluctuations is 
restricted by the increase in the free energy (areaxcr) of the domains.^ The resulting 
probability distribution is a "smeared" version of the thin-wall case, with the pure 
phase 5-peaks rounded into Gaussian distributions. This is shown as thick lines in 
Figs. [I and |. 

^In order to be able to properly define the domain shape in fluctuating, real systems, one should 
coarse-grain the order parameter up to, at least, the length scale given by the bulk correlation length. 
In the calculations in this paper we use only volume averages of the order parameter, and the coarse- 
graining is of no consequence. 
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The probabihty distribution P{6) is to be measured with Monte Carlo simulations. 
Due to the enormous suppression of the critical droplets, standard Monte Carlo sim- 
ulations, where configurations are sampled with weight p oc exp[—H], are completely 



inadequate for the task. However, the multicanonical Monte Carlo method |2T[ has 
been developed for just this kind of task; this will discussed in Sec. ^I2| . From the 
measured probability distribution P{9) oc exp(— F(6')) at small supercooling we 
can measure the critical droplet free energy 

Fc( W) ^ F{ec) - F(^symm.) = " log f/^^ , (2.5) 

where 9c is the local maximum of the free energy, and ^symm. is the location of the 
symmetric phase minimum. Since one should remain a safe distance away from the 
cylinder and slab regions, X should remain well below ~ 0.15. This limits the range 
of 5rn?' which can be probed at a given L, or sets what L is required for a given dm^. 
We note that Eq. ( p.5| ) resembles the widely used histogram method for calculating the 
planar interface surface tension a, proposed by Binder [^. In some cases Eq. ( ^.51 ) 
can already give a fair estimate of the droplet nucleation rate: 

r oc [time]"^ [length]"^ exp(-Fc;) , (2.6) 

where [time] and [length] are suitably chosen dimensionful constants. Unfortunately, 
Eq. ( p.6|) tells us nothing about the constant of proportionality, and the dynamical 
information is missing completely. The result also explicitly depends on our choice of 
the order parameter. Thus, the rate can be off by orders of magnitude. Furthermore, 
as opposed to Binder's method, which gives the surface tension a exactly as ^ cxd, 
Eq. ( |2.5| ) does not have a well-behaved infinite volume limit with fixed 5m?. The main 
problem is that, for a fixed supercooling and droplet size, we have 9c— Osymm oc l/V, but 
the width of the Gaussian-like symmetric phase peak is oc Thus, the symmetric 

phase fluctuations will "swallow" 9c when the volume is too large! In practice, this 



does not appear to be a problem for suitably chosen order parameters (see Sec. 4.1 



2.2 Rate of droplet nucleation 

Although the droplet nucleation rate is inherently a non-equilibrium problem, in a finite 
(but large) volume it can be formulated fully in thermal equilibrium. Let us take an 
ensemble of configurations in thermal equilibrium at small supercooling (5m^, and let 
each of the configurations evolve according to some dynamical prescription. Naturally, 
we require that the evolution preserves the equilibrium ensemble, and we consider here 
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only Markovian dynamics which satisfy detailed balance 



T[{(j)b, -Tib) (0a, -VTa)] 



exp[i7(0b,7rb) -i7(</.,,7r,)]. (2.7) 



Here are field variables, vr the momentum variables (which may not exist in the 
dynamics at all). if(0, vr) = if(0, —it) is the energy of the state (scaled by T to be 
dimensionless) , and T(a b) is the transition probability (along some trajectory) 
from state a to b. All standard dynamics we know of obey this condition (for example, 
Hamiltonian, Langevin, Glauber, heat bath). 

We specialize to supercooling,]^ so < m^, and at any small time interval the vast 
majority of the configurations remain in the broken phase, a much smaller fraction 
are in the symmetric phase, and a tiny number of configurations are in the process of 
evolving from one phase into another. In equilibrium the rate must be equal in both 
directions. Let us choose two order parameter values 9s and 6*^ on the symmetric and 
broken phase sides of 9c, chosen so that the constrained free energy F{9s), F{9b) -C 
F{9c)', see Fig. 0. If we now take a configuration from the thermal distribution with 
order parameter restricted to 9s and let it evolve in time, due to the steep slope in 
the free energy, it is almost certain not to evolve to 9c in the "medium term," i.e. in 
timescales shorter than the very long tunneling time l/TV. Similarly, a configuration 
at 9b is overwhelmingly more likely to fall down into the broken phase than to go to 
the critical droplet. 

During the tunneling the configuration must evolve along a continuous trajectory 
through the order parameter values 9s ^ 9c —* 9b- We shall focus on the piece of 
the real time evolution trajectory between the last time the order parameter had the 
value 9s, before reaching 9b for the first time. Since the time evolution of the order 
parameter may be "noisy," this trajectory may cross the value 9c several times in a 
short time interval. We note here that the values of 9s and 9b are not critical, as 
long as the free energy condition F{9s,b) ^ F{9c) is satisfied. For example, we could 
choose ^5 to be well within the bulk of the symmetric phase; the point is that the 
trajectory 9s ^ 9b can be used to identify the tunneling events. 

We define the symmetric phase probability Psymm. to be the integral of the canonical 
distribution P{9) oc log[— ^(6^)] up to up to the critical droplet value 9c, 

Psymm. = / d9P{9) , (2.8) 



with 6'mm the smallest possible value of 9. The choice of 9c here is for convenience; 
since by far the dominant contribution to the integral comes from the central peak of 
^The case of superheating, tunneling out of the broken phase, is a trivial extension. 
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the symmetric phase, choosing some other reasonable value in the neighborhood of 6c 
would give only exponentially suppressed corrections. 

Let ns{t) be the probability that a configuration which was in the symmetric phase 
at time t = still remains there at time t, normalized so that ns{t = 0) = Ps- The 
tunneling rate F is defined through 

TVn,{t) = ^ ns{t) = risiO) exp [-TVt] . (2.9) 

This immediately suggests that we can measure F by simply measuring the average 
time r it takes for the symmetric phase configurations to tunnel: F = 1/ (Vt). However, 
when the rate is very small, as it is in our case, this approach is utterly impractical.|^ 
Let us now consider times which are much shorter than the tunneling time, t <C 
l/{Vr) (this will not restrict the validity of the final results). In this case the symmetric 
phase consists almost entirely of configurations which were there already at t = 0. We 
can substitute ns{t) ^ ns{0) = Pg on the l.h.s. of Eq. (|2.9| ), and dns{t)/dt becomes 
the time-independent fiux of configurations from the symmetric to the broken phase. 
The central idea of the method described here is the realization that, since all of the 
tunneling trajectories must go through 9c-, we can calculate the fiux dng/ dt by sampling 
configurations and trajectories directly at ^ = 9c- Labelling the configurations with 
symbol and the trajectories through with a^, the fiux can be written as 

^ = I d<pda^ Ae-^(^'"^) 5{9{<P) - 9c) x x (2.10) 

Here Aexp — /f(0, a^) gives the thermal equilibrium probability of a configuration 
and a trajectory which goes through it. In other words, we want to average over 
all configurations such that ^(0) = 9c-i and all trajectories which go through the 
configuration 0, with the correct, thermal weight. A is a normalization constant, 
which satisfies 

" rf0rfa<^Ae-^('^'°^)5(e(0) -9)= P{9), (2.11) 



where P{9) is the order parameter probability distribution. 

The dynamical information is contained in the last two terms of the integrand in 
Eq. (|2.10| ). The derivative d9{a(i,)/dt is evaluated along the trajectory at the point 
where crosses the pivot configuration 0. Its purpose is to cancel the Jacobian factor 



Strictly speaking, at extremely long times Eq. ( 2.9 ) is not correct; it neglects the multiple tunneling 
trajectories symmetric — > broken symmetric etc. Indeed, in the limit t —> oo, n{t) « ^'symm • 
However, since the rate of the the inverse process is much smaller than VT, this becomes significant 
only at times which are much longer than the tunneling time (TV)^^. 
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arising from the (5-function: along the trajectory, 



5(9 - 9, 



c) 



d9 ^ 



dt 



S{t-tc), (2.12) 



where tc is the time when the configuration crosses the point 9c- Thus, it turns the 
integrand from a probabihty density of configurations at 9c to a flux of trajectories 
through the value 9c- This can be understood in another way: if d9/dt is small for 
trajectory a, a conflguration evolving along it loiters near 9c longer than a conflguration 
on a trajectory where d9/dt is large. Thus, a "slow" trajectory contributes more to 
the probabihty density e~^5{9 — 9c) than a "fast" one; and the inclusion of d9{a)/dt 
exactly cancels this bias. 

Finally, in order to count only trajectories which really tunnel from the symmetric 
phase to the broken phase, we deflne the projection Sg_,Q = 1, if a tunnels from 9s to 
^B'l ^s^B — otherwise. As discussed above, we deflne the tunneling trajectories to 
be ones which evolve from 9s through the conflguration (p to 9b- The result is that the 
integral in Eq. ( p.lO| ) evaluates the number density of tunneling trajectories through 
9c, which is exactly the rate we are after. 

During the evolution from 9s to 9b the tunneling trajectories may cross 9c several 
times (always an odd number), and these trajectories are counted by the integral (|2.10D 



at every crossing. There is no overcounting, however: if there are (2n+l) crossings, in 
(n+1) of these d9/dt is positive, and in (n) it is negative, corresponding to crossings 
left— >-right and right— i>left, respectively. The absolute magnitude of the integrand is 
equal at every crossing of 0avC along the trajectory: it follows from Eqs. ( |2.71 ), ( p.lO|) 



and Eq. ( p.l2|) that the probability of choosing a crossing conflguration (pa and a tra- 



jectory which evolves to another crossing conflguration (pa — *■ (pb, is equal to choosing 
one evolving (pb — >■ (pa- Thus, n positive and negative contributions cancel, and the 
trajectory is counted only once.[^ 

This also implies that the dynamical factors in Eq. ( p^.lO[ ) can be rewritten in the 



following form, which is much more suitable for numerical work 

d9{a 



de{a) 1 
~ir ^'-^ 2 



dt 



X . (2.13) 

crossings 



Here N^^^^s,ing is the number of times the trajectory a crosses 9c, and ^J^nnei = 1 
trajectories which tunnel either 9s — > 9b or 9b — >■ 9s'- since the whole ensemble is in 
thermal equilibrium, the flux of trajectories from the symmetric to the broken phase 
must equal the inverse one. We can take an average over the two sets of trajectories. 
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See also the discussion in the appendix A of Ref. [Qi 
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hence the factor 1/2 in ( ^^.131 ). Using the expression ( |2.13| ) makes the integrand in 



Eq. ( |2.1CI| ) positive definite, which improves the convergence of the numerical Monte 
Carlo integration. 

In Eq. ( p. 10 ) we implicitly assumed that the trajectory is differentiable in t. This is 



true for Hamiltonian type dynamics, but not for "noisy" time evolution, like Langevin 
or Glauber dynamics. In these cases we must substitute the derivative with a finite 
difference: 

and the trajectory must be sampled with the same frequency At when the crossings 



in Eq. (|2.13|) are counted. Due to the noisiness of 0{t), the smaller At is, the more 
crossings of 9c are potentially resolved. However, this is compensated by the increase in 
A9/ At. For example, in Langevin dynamics, when At is small enough, the trajectory 
9{t) behaves almost like a Brownian random walk in the proximity of 9c- In this case 
(|A6'/At|) ~ (At)-i/^ and (iVcrossings) ~ {At)-^'^. Thus, neither of these terms have a 
At ^ limit, but the product in Eq. (|2.13|) has. 

The use of finite time differences arises automatically in numerical time evolution, 
also in Hamiltonian evolution, and we shall use this formulation from now on. It 
is natural to use the same At in sampling the trajectory as is the step size of the 
time evolution; however, this is not necessary and the sampling interval can be longer. 
However, it should not be so long that the free energy F{9) changes appreciably in one 
interval. In practice this is not a problem. 

2.3 Monte Carlo evaluation 



In practical Monte Carlo simulations we evaluate Eq. ( p. 10 ) as follows: let us select 



configurations from the thermal ensemble, p oc e~^, but with the order parameter 
restricted to a narrow range \9 — 9c\ < e/2. We can now construct trajectories (one or 
more for each configuration) going through by evolving the configuration both forward 
and backwards in time, until both the forward and the backwards ends of the trajectory 
reach order parameter values 9$ or 9b- The trajectory contributes to tunneling only if 
the ends of the trajectory are on different sides of 9c- Thus, the Monte Carlo version 
of Eq. ( p.lO| ) becomes 

■ A9ia) 



TV = 

where we have defined = ^runnei/^Sossings. 




xd"), (2.15) 



Pc = -J5 / . (2-16) 

^-^ symm. ''Sc~2 
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Here we have substituted the 5{6 — ^c)-function in Eq. (|2.1CI|) with a small but finite 



width e; in practice, it is much easier to generate configurations in this ensemble 
than with a sharp 5-function. The derivative d9/dt is again evaluated at the initial 
configuration 0, and the expectation value (■) is over the configurations and trajectories. 

The evolution backwards in time means that we invert the (possible) initial momenta 
TT — >• — TT and evolve the system using the standard equations of motion; the time is just 
interpreted as —t. If the dynamics does not involve momenta (for example Langevin 
or Glauber), then the evolution backwards is just another realization of the standard 
forward evolution. The detailed balance condition ( |2.7| ) guarantees that the trajectory 
generated by this backwards evolution can be re-interpreted as a representative forward 
evolving trajectory. 

The expression (p^.lSD is straightforward to use in simulations. However, we simplify 



it here a bit further: for the physical process we are considering in this work, it turns 
out that |A^^(a;)/At| is almost completely uncorrelated with the "global" properties 
of the trajectory a, that is, how many crossings of 9c it makes, or whether it finally 
contributes to tunneling. This implies that to good approximation we can evaluate the 
expectation values of the two terms inside the angle brackets separately: 



TV 




A9 



At 



(d) . (2.17) 



Both of the expectation values above are thermal averages over all trajectories which 
go through 9c- This is the form of the rate equation given in Ref. [0. The advantage 
here is that (|A6'/At|) is a 'local' quantity, and its evaluation does not require the full 
knowledge of the trajectory. Indeed, for the cubic anisotropy model, with the order 
parameter and the real-time dynamics we shall consider in this work, it is possible 
to evaluate it analytically, see Sec. The decomposition of the expectation value 
in ( [^.17D is a natural consequence of the fact that the critical droplets are large ob- 
jects, and thus must involve a large number of microscopic variables. The evolution 
of the UV-modes of these variables is largely uncorrelated across the droplet. The 
decomposition holds strictly in the zero lattice spacing limit, where \d(f)l^/dt\ is UV 
dominated. However, if one applies the method described here to tunneling processes 
where only relatively few microscopic degrees of freedom participate, one cannot rely 
on the decomposition and the equation ( |2.15| ) must be used instead. 



2.4 Independence on the order parameter 

In the discussion above we made extensive use of the critical droplet order parameter 
value 9c, where the free energy F{9) has a local maximum. However, we could have 
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used any other value of 6 in the neighborhood of the free energy maximum — we 
only have to make sure that the whole of the bulk of the symmetric phase probability 
remains on the symmetric phase side of it. Intuitively this is clear. What Eq. ( p^.lOD 
actually measures is the number of trajectories going from 6^5' to Ob] and this must 
be the same, no matter what point in between we decide to intercept them. For 
example, we may choose a new value 6' , somewhat to the symmetric phase side of 6'c, 
to be our new starting point. The probability factors in the integrals ( 2.10 ) or ( |2.15 ) 



increase (the free energy is lower), but this is exactly compensated for by the decrease 
in the average value of the projection 5'^^^^^^ most of the trajectories which intersect 
9' will both start and end in the symmetric phase and hence do not contribute to 
the tunneling. However, if we want to get as accurate a result as possible for a given 
amount of cpu-time, it pays to choose the starting value as close to 6c as possible. 

Very significantly, the rate calculated with (p.l0|) -( p.l7|) is independent of the choice 
of the order parameter itself. This is in contrast to the probability distribution P{0) 
alone (which can be used to calculate the free energy of the critical droplet, see Sec. |2?T| ). 
The requirement the order parameter has to fulfill is that it separates the symmetric 
and broken phases and resolves the potential barrier between them to a sufficient 
degree, in order for one to be able to choose the points corresponding to 6s, 9c and 
9b- Provided that the hierarchy in free energies F{9c) ^ F{9s),F{9b) is satisfied, 
the order parameters yield equal results, up to exponentially suppressed corrections, 
of order exp[(F(^5) - F{9c))] < 1. 

However, not all order parameters are created equal: the best parameter has as small 
fluctuations in the bulk phases as possible, while separating the two phases as strongly 
as possible. Large fluctuations will reduce the resolving power of the order parameter: 
for example, we may have a droplet which is slightly smaller than the true critical 
droplet, but a fluctuation in the bulk may still bring the order parameter up to 9c- 
These fluctuations will typically increase the probability P{9c), and again, this increase 
is exactly cancelled by a decrease in 5^nnei- ^^e trajectories evaluated from these con- 
figurations are less likely to lead to tunneling. Thus, the better the order parameter, 
the larger the fraction of trajectories which lead to tunneling. Especially important is 
that the fluctuations of the order parameter in the symmetric phase should be very 
small, since more than 85% of the lattice volume in a critical droplet configuration will 
remain in the symmetric phase, outside the droplet. 
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3 Cubic anisotropy model 



For concrete calculations in this work we use the cubic anisotropy model in 3 dimen- 
sions. The reasons for choosing this are the following: (1) it has a radiatively generated 
first order phase transition, the strength of which can be adjusted continuously; (2) it 
is a (superficially) simple model and relatively easy to study on the lattice; (3) it is a 
continuum field theory, and we can test whether our nucleation calculation gives con- 
sistent results as the continuum limit is taken on the lattice. Furthermore, recently the 
validity of hanger's theory of nucleation has been questioned in radiatively generated 
transitions, using exactly this model as an example |T^ . 



3.1 Definition of the model 

The cubic anisotropy model is a model with two scalar fields (pi and (j)2, with quartic 
interactions respecting the discrete symmetries (pi ^ 02 and (pi ^ —(pi- The ther- 
modynamics of the model is defined by the partition function (We use here "particle 
physics" units where [length] = [time] = [energy]"^ = [temperature]""^). 



V(j)iV(P2exp-Ho[(P], (3.1) 



where the 3-dimensional action is 



HM] = / d^x 



a=l,2 \ 1=1 




(3.2) 



The scalar fields (pi, (p2 have dimensions of (length) ' , and the coupling constants A ~ 
(length)"^. In the case that the 3 dimensional theory arises from dimensional reduction 
of a 3+1 dimensional, relativistic scalar field theory, the 3-D quantities are related to 
the 3-t-l dimensional quantities through Hq = H^d/T, (p = 04dT^/^, A = A4dT, and 
= miY) + 0{T'^). Besides rescaling the 3 dimensional parameters, changing the 
temperature of the 3+1 dimensional system adjusts the value of m^, which is why we 
often equate lowering with cooling. 

The mass parameter drives the transition: when has a large enough positive 
value, the system is in the symmetric phase {(pi) = {(p2) = 0, whereas with large and 
negative we are in the broken phase where \{(pi) \ + \ {<p2)\ > 0. At the mean field 
level ("tree level" in particle physics parlance) the transition between the symmetric 
and broken phases is of second order and happens at = 0. Also at the mean field 
level, the stability of the theory requires Ai > and A2 > — Ai/3. We can immediately 
recognize some special values for the couplings (Ai, A2): first, when A2 = Ai/3, the 4th 
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order term in the energy functional becomes (Ai/24)(</)^ + (j)^'^. Thus, the theory has 
then an 0(2) symmetry under rotations of the fields (0i, 02), and as we vary m^, there 
is a second order 0(2) phase transition. Another special point is A2 = 0, where we 
have two uncoupled 0^ models. Now the transition is again of second order, of Ising 
model universality. We also note that the redefinition of the fields 

, , , . 1 
72 



1,</^2J —1</^1 + 02,01 - 02) (3.3) 



maps Eq. (|3l^) onto itself, with 



m^^m\ Ai^^(Ai + 3A2), A2 ^ ^(Ai - A2) • (3.4) 

The 0(2) invariant point A2 = Ai/3 is a fixed point of this transformation, and the 
region A2 > Ai/3 is mapped on A2 < Ai/3. This implies that the point Ai = A2 is 
equivalent to A2 = 0, and it too corresponds to two uncoupled 0^ theories with Ising 
universality. 

The behavior of the system at general values of the couplings has been studied by 



renormalization group methods by Rudnick [jTI| and Arnold and Yaffe [|T^]. Their 
analysis indicates the following behavior: 

1. If < A2 < Ai, the theory will flow to the 0(2) fixed point at large distances, and 
so will have a second order transition with 0(2) criticality. In the broken phase 
there is a Goldstone mode, corresponding to rotations of the (0i, 02)-vector. 

2. If A2 = or A2 = Ai, we have two uncoupled 0^ theories, and the transition is of 
Ising type. 

3. If A2 < or A2 > Ai, the theory does not flow to any weakly coupled infrared 
stable fixed point, and so might be expected to have a first order phase transition. 
This is indeed the shown by perturbative analysis O, O] and lattice 



simulations 13 



The first order transition in the last case is generated by radiative effects; in order 
to observe it perturbatively, one must consider the first loop correction (Coleman- 
Weinberg effect |2^ ) . Using the e expansion around d = 4 — e dimensions, Arnold and 



Yaffe |T2[ calculated the effective potential up to next-to-next-to-leading order in e. In 
the limit A2 ^ Ai, the transition becomes strong and the perturbative analysis should 
be accurate. There are 4 degenerate broken phases in the region of the first order phase 
transition. If we choose A2 > Ai, in the broken phase the (0i,02) -doublet acquires an 
expectation value along one of the principal axes: (±t>, 0) or (0, ±f ). 

In this work we are interested in strong first order transitions. Thus, we choose 
A2 = 8A1, which is seen to give a very strong transition. 
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3.2 Cubic anisotropy model on the lattice 



We discretize the theory on the lattice using an C(a^) accurate lattice action, where a 
is the lattice spacing. The action can be written as 



Hn 



E 



(-0iV 



L01 



L(P2 



(m 



5m? 



? + 02 + 
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(0^ + 0^ + 



(A2 + (5A2) 

4 



(3.5) 



where we have rescaled the lattice fields and couplings to be dimensionless: (/)iatt. = 
0cont.a^^^5 Aiatt. = Acont.a, and mf^^^ = m^^^^o?. 5-terms and Z-factors are additive and 
multiplicative renormalization constants, respectively. They are necessary to match 
the behavior of the lattice theory to the continuum theory without O(a^) errors, and 
their calculation is detailed in the appendix. Eliminating such errors also requires a 
lattice Laplacian which includes next-neighbor couplings. 



1 5 4 
VMx) = -y0(x) + -E 



x + i) + (h[x — i) — 



{x + 2i) + d)ix-2i 



(3.6) 



Furthermore, the operator insertion {(fp') = {cfl + (f)\), which we will need in what 
follows, is additively and multiplicatively corrected. 



Ct(0^)contin 



latt 



(3.7) 



3.3 Real time dynamics 

To define the nucleation rate, we must specify a dynamical prescription for the time 
evolution. The rate will depend on the choice, which, after all, determines what physical 
system we are describing. Any prescription which preserves the canonical distribution 
is permitted. The most common evolution prescriptions in field theories are the Hamil- 
tonian and the Langevin dynamics. Both are perfectly valid evolution dynamics for 
classical thermal field theories; for our case the Hamiltonian dynamics has the further 
conceptual advantage that it also describes the evolution of a quantum theory at high 



temperatures to leading order in coupling constants [^, ^ . 

In this work we shall consider a one-parameter class of evolution dynamics, the 



Hamiltonian stochastic equations [^6 



dt(pa{x,t) 
dtna{x,t) 



1Ta{x,t) 



+ -fTTa{x,t) +^a{x,t) 



(3.8) 
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Here Hq is given by Eq. ( |3.2| ), tTq are the momentum variables of fields 7 > is 
a "friction coefficient," and is Gaussian noise, which satisfies the stochastic 

condition 

iU^, tM^', t')) = 2^6at6''ix - x')5{t - t') . (3.9) 

The equations of motion ( p.8| ) thermalize the system (if 7 > 0), that is, at large times 
the probability distribution approaches the Gibbs distribution p oc exp(— vr]), 
where 

i/[0,7r] = ld'x-inf + nl) + So[(l)]. (3.10) 

Thus, the evolution also preserves the correct thermodynamics of the "static" theory, 
Eq. dO). 



If we choose 7 = in Eqs. (3^) the evolution becomes the standard canonical 
Hamiltonian evolution, with canonical momenta tTq and conserved energy H. On the 
other hand, if 7 is very large, the equations of motion can be reduced to the familiar 
fully diffusive Langevin equation for 0, simply by eliminating vr and neglecting the 
term (small in this limit) d^cp. Thus, by adjusting 7 we can continuously adjust the 
coupling of the system to an external "heat bath." 

For the case of droplet nucleation at finite volume the presence of noise improves the 
finite volume and lattice spacing scaling of the dynamics. This is because the transition 
has a large latent heat. A growing/shrinking droplet will release/absorb a significant 
amount of energy. Under microcanonical Hamiltonian dynamics the temperature in a 
finite system will correspondingly increase/decrease. In an infinite volume, the released 
energy would be transported rapidly away from the proximity of the droplet by hydro- 
dynamical fiows0, and the temperature would remain close to constant. As the lattice 
spacing becomes smaller, the heat capacity of the system grows, and the latent heat 
again becomes less important. But at finite lattice spacing and volume, the latent heat 
remains indefinitely under Hamiltonian dynamics, modifying the time development. A 
finite 7 in Eqs. ( p.8|) absorbs the latent heat, allowing the system to more rapidly ap- 
proach the infinite volume behavior. It does this efficiently for 7 > 1/L, in which case 
the momenta thermalize completely in time r <L. In Sec. ^ we check the dependence 
of the nucleation rate on 7, and we find that we can vary 7 over a large range without 
significantly affecting the results. This means that, in practice, the nucleation rate has 
a fairly weak dependence on the specifics of the real time dynamics. 



^^Hcat flow behaves very differently in a plasma without nonzero conserved quantities than in more 
familiar condensed matter or convective settings, see for instance p^. 
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3.4 Real time evolution on the lattice 



The discretization of the equations of motion (|3.8| ) requires some care, because certain 
errors lead to time evolution which respects a slightly different ensemble of configura- 
tions than the Gibbs one. This is important, because slight thermodynamic changes 
dramatically affect the critical bubble behavior. The easiest safe update at nonzero 7 is 
to apply Hamiltonian evolution, stopping every At ^ a to apply a partial momentum 
refresh, 

7r,(x, t + 0) = (1 - 7r,(a;, t - 0) + er]a{x, t) , (3.11) 

with = (1 — exp(— 27At)). This refresh exactly preserves the Gibbs distribution. We 
perform At of Hamiltonian evolution by a single step of the fourth order Runge-Kutta 
algorithm ||29], which has errors suppressed by {At/a)^. This procedure reproduces 
Eq. ( |3.8| ) with errors of order 0([At/a]2), 0{e^), but what is more important, it main- 
tains the Gibbs distribution to 0([At/a]^) errors. We use At /a = 0.05, which proves 
more than adequate, unless we want to study extremely large 7. 



4 Numerics and thermodynamic results 
4.1 Order parameter 

The cubic anisotropy model has a true order parameter, which measures the breaking 
of the global sjTiimetry: 

0av = \/(^^j^T(^, ^l = ^,J2h- (4.1) 

Here A^^ is the volume in lattice units. However, a better observable for our purposes 

0L = ^E(0? + 02). (4.2) 

This is clearly not an order parameter in the sense that it would be zero in one phase 
and non-zero in another. However, it has different expectation values in the symmetric 
and broken phases, which is sufficient for our purposes. The reasons why we prefer 
over 0av are that (a) it is much more economical to use, and (b) it resolves the critical 
droplet more accurately. 

Let us first consider the cost. Note that (pl^ is multiplied by in the Hamiltonian 
( |3.2| ). Assuming that we have measured the order parameter distribution at = ml, 

Pm0l) oc I rf0'exp[-if^.(0')] 5i<l>Z - 0L) , (4.3) 
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we obtain the distribution at any other simply by reweighting: 



(4.4) 



The reweighting can be done without any significant loss of accuracy: if we have 
determined P^2 to a good accuracy in some ^^^-range of interest, reweighting gives P^2 
with the same relative error SP{(f)l^)/P{(j)1-^) in this range.Q Since the calculation of 
-P(0av) comprises the major part of the total computational effort, it is quite significant 
that we can perform the measurement once, and use it for all values. Of course, 
separate calculations are still required for different lattice volumes and lattice spacings. 

As discussed in Sec. ( p.4[) , in order to be able to resolve the droplet as well as possible, 
the fluctuations of the order parameter in the symmetric phas^ should be as small 
as possible — this is because more than 85% of the volume of the lattice is in the 
symmetric phase in critical droplet configurations. Here (pl^ is a clear winner over 0av 
This can be seen by comparing the (width)^ of the symmetric phase distributions, i.e. 
the susceptibilities: 



Here is the (mass)^ of in the symmetric phase. Since the tree-level transition 
happens at = 0, the value of at the transition must be small, of order A^. The 
relevant quantity at the transition is the width of the symmetric phase distribution 
compared against the broken phase expectation value; using 1-loop effective potential. 



where v = (0av) in the broken phase. Thus, the symmetric phase fluctuations of 0av are 
more than 10 times larger than 0^^, when normalized to the broken phase expectation 

^■^We can also reweight P(0av) with respect to m? (or any other order parameter distribution, for 
that matter); however, in this case the errors increase dramatically when is outside a narrow range 
around the value where the original measurement was done. 

^"^In general the argument applies to the phase the system is tunneling out of. In particular, 
would be a very bad order parameter if we were studying tunneling from the superheated, broken phase 
into the symmetric one; we would want to choose instead some measurable with small fluctuations in 
the broken phase. 




(4.5) 



(4.6) 





(4.7) 
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values. Naturally, the final word comes from the lattice simulations, where (psx is indeed 
observed to give a ~ 10 times wider symmetric phase peak than 0^^. The width of the 
symmetric peak of the the 0av-distribution is actually larger than the critical droplet 
value 0av,c for the parameters we are using in this work. Thus, 0av,c does not have the 
necessary resolution power, while we can use without problems. A higher power 
of 0, for example, (f)%j, might give an even narrower symmetric phase peak on a coarse 
lattice spacing, but in the continuum it would have a UV divergent susceptibility, unlike 
so taking the continuum limit might prove problematic. Such order parameters 
also fail the economy criterion. 



4.2 Multicanonical method 

In order to be able to calculate the droplet nucleation rate, our task now is to de- 
termine P(0^y) in the droplet regime, X < 0.15. Since the probability density can 
vary by ~ exp(lOO) or more in this range, normal Monte Carlo simulations, where the 
configurations are sampled with probability p oc exp —H, are completely out of the 
question. The multicanonical method was originally developed for just this kind 
of problems. In multicanonical sampling the configurations are chosen with a modified 
weight 

p^ucacxexp[-if„2 + iy(02j]. (4.8) 

is a weight function, which is carefully chosen so that the probability distri- 
bution of measured from the multicanonical Monte Carlo simulation, Pmuca(0av)) 
is approximately constant in the range of interest. The canonical distribution is then 
obtained by reweighting with the weight function: 

0(<^L)oce-^(*^v)p^_(<^2^). (4.9) 



This can be further reweighted to other values of using Eq. ([4.4|) . 

An optimal choice for W is — logPm2(0av)' which is just the quantity we want to 
calculate with the multicanonical simulation! The weight function need not be ex- 
actly the ideal one, but it must not deviate from it too much or the multicanonical 
simulation becomes very inefficient. This chicken-and-egg problem can be resolved by 
first calculating W with an automatic iterative feedback procedure PO, PTI. We use a 



variation of the one presented in Ref . |^ . This method yields progressively better ap- 
proximations for W , until sufficient accuracy is reached. The final W can then be used 
in a multicanonical simulation, which finally gives us -P(0av)- ^ detailed description of 
the application of the multicanonical method to the problem of the bubble nucleation 



in the electro weak theory can be found in |jT9| 
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In the multicanonical phase of the simulations we are at hberty to choose as efficient 
an update algorithm as possible. (In the real-time runs the evolution equations fix 
the update.) We use a mixture of site by site heat bath and a non-local multi-grid 
over- relaxation. 



4.3 Results: thermodynamics 

Before tackling the droplet nucleation rate calculation, we want to study basic thermo- 
dynamic quantities associated with the phase transition: the transition "temperature" 
= m^, the "latent heat" ^, and the tension of the interface between the symmet- 
ric and broken phases a. We want the above quantities in the true thermodynamic 
limit, that is, we extrapolate to infinite volume and to zero lattice spacing. All of the 
simulations have been done with A1/A2 = 1/8, which guarantees a strong first order 
transition. 

In this part of the analysis we use up to 6 different lattice spacings a\2 = 1.5, 2, 2.5, 
3, 3.5 and 4 (all dimensionful quantities are given in terms of A2). For every lattice 
spacing we perform simulations using a series of lattice sizes; all in all, we have 64 
simulations at different lattice spacings and volumes. Since the transition is strong, all 



of the simulations described here are multicanonical (see subsec. [4. 21) , with the weight 
function optimized for the whole order parameter range from the symmetric to the 
broken phase. 



The transition point m^: We define the transition value of to be the value where 
the two peaks of the probability distribution P(0^y) have equal weight, i.e. the value 
where the volumes of the symmetric and broken phase peaks are equal (see Fig. |^). 
More precisely, the symmetric and broken phase probabilities are the integrals of the 
distribution to the right and to the left of a separating value between the peaks; since 
the probability is exponentially suppressed here in sufficiently large volumes, the precise 
choice of the separation value is exponentially unimportant. For each lattice spacing 
a, we determine the transition point using a series of volumes and extrapolate the 
results to infinite volume. The extrapolation is linear in 1/V, as long as the volumes 
have similar geometry; an example of this at a = 3/A2 is shown on the left panel of 
Fig. I 

The infinite volume points are in turn extrapolated to the continuum limit a — > 
0. This is shown on the right panel of Fig. ^. Since we know the additive mass 
counterterms only to order 0(a°) (see appendix we extrapolate to the continuum 
limit using an ansatz cq + cix + C2X^ + c^x^, where q are fitted constants and x = a\2. 



22 



10" 



X,a=2.5, y=22. 



10" 



-Q -20 

o 



10" 



10" 




Figure 3: Reweighting of the probability distribution P{(j)1^) to different values of m^. 
The transition value is the one where the symmetric and broken phase peaks have 
equal volume. 



The final result is 

mlifi = \2)/\l = 0.0096 ± 0.00022 (4.10) 

We note that the supercooling Sm^ can be determined much more accurately than m^: 
in differences of the additive counterterms cancel, and the leading errors behave as 
O(a^). At no stage of the nucleation analysis do we need to know the absolute value 
of m^. 

In the fit the coefficients Ci and C2 are actually consistent with zero, suggesting that 
the unknown counterterms are small. We also show a fit with these coefficients set to 
zero. 



The surface tension cr: The surface tension is the free energy per unit area of 
interface. As discussed in Sec. |2.1| , when = the suppression of the probability 
in the mixed phase (in large enough volumes) is caused by the existence of phase 
interfaces. Approximately half-way between the symmetric and broken phases, the 
mixed phase configurations consist of slabs of symmetric and broken phases, separated 
by two approximately planar interfaces. Thus, the interface tension a can be obtained 
from m 

(7 = — lim — — — log — . (4-11) 

V^oo 2LxLy -Ppcak 
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Figure 4: Transition value of extrapolated to infinite volume at = 2 (left), and 
to zero lattice spacing (right). The solid line shows a 3rd order polynomial fit, and the 
dotted line is a fit of form cq + Cs{\2aY. 
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Figure 5: Surface tension (t/X^ extrapolated to infinite volume at = 3 (left), and 
to zero lattice spacing (right). 

Here we have assumed that the interfaces are oriented parallel to (x, y)-p\ane. As for 
other quantities we have scaled T out of a; a more conventional definition is cTc = crT, 
so that (Tc X area is the free energy of the interface (rather than the scaled free energy we 
use in this paper). At the "equal weight" point the symmetric and broken phase heights 
are not usually equal; in this case Ppeak is the average of the two peaks. However, here 
we determine a from distributions which have been reweighted to "equal height" m^; 
in our case, the two definitions give completely consistent results. 

The extrapolation to infinite volume in Eq. ( [4.11| ) can be substantially improved by 
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using finite volume corrections 
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aa^ = ^ log (^) logn. - \ logn. + \g + const.) . (4.12) 

z^iOx'i'y -I slab '''x''"y ^ ^ ^ 

Here Li = riia and ai = oa^ is the tension in lattice units. The function G interpolates 
between lattice geometries. For cubical volumes G = log 3, while for long cylinders 
G = 0. In our analysis we use strongly cylindrical boxes, ^ = Ly. This 
guarantees that the interfaces form along the (x, y) plane, and it also separates the 
two interfaces much farther than a cube of the same volume. On the left panel of 
Fig. ^ we show the extrapolation of cr/Ag (corrected with Eq. (|4.12|) ) to infinite volume 
at lattice spacing = 3. The continuum limit extrapolation is shown on the right 
panel. Because of our lattice improvement the leading errors are O(a^). We obtain a 
good fit using a third order ansatz a{a) = cr + c^a^, with the result 

(t/A^ = 0.0050 ± 0.0002 . (4.13) 



The latent heat i: Since is our temperature parameter, we define the latent 
heat using 

^) = ^((Obroken ' (Osymm.) ■ (4.14) 



The difference between the symmetric and broken phase expectation values is readily 
measurable from the distributions like the one in Fig. ^. Our extrapolation procedure 
is the same as for a: 

£/A2 = 0.243 ±0.004. (4.15) 



5 Results: the nucleation rate 

We evaluate the rate of the droplet nucleation using Eq. ( |2.17| ). To recap the discussion 
in Sec. 0, the calculation consists of two separate stages: 

(1) the measurement of the probability distribution of the order parameter, P((j)l^), 
using multicanonical simulations, and 

(2) the real-time evolution of the critical droplets, which gives us (|A0^y/At|) and (d). 
The lattice spacings and sizes used are shown in Table |l[ For each of the lattices we 

calculate -P(0av) once; as explained in Sec. (|4.1| ) and (|4.2|) , the original distribution can 
be reweighted to different values of supercooling Sm^. On the other hand, the real-time 
trajectories have to be calculated separately for each supercooling. 
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a\2 


size / 


evolution at 5m^/A| 


4.0 


60^ 


0.001156 


3.5 


703 




3.0 




0.0007-0.00125 


2.5 


1003 


0.00094 


2.0 


1203 





Table 1: The lattice spacing and size where the droplet nucleation rate has been 
calculated. For each lattice the probability distribution P(0av) been calculated 
once with a multicanonical simulation. The real-time evolution of the critical droplets 
are calculated using the supercoolings shown on the third column. 



5.1 The probability of the critical droplets 

In order to fit the droplets comfortably inside the lattice, the lattice volumes have to 
be substantially larger than the volumes used in Sec. ([4.3|). However, as discussed in 
Sec. (|2.1|) , in this case we have to calculate the probability distribution P{(p1^) only in 
the "droplet branch" of the distribution, < [0^v - 0av,s] < 0.15 x [0^^ ^^ - (f)l^J. This 
guarantees that the distance of the droplet from its periodic copies is large and we 
are a safe distance away from the "cylinder" branch of the distribution. Furthermore, 
the restricted range of reduces the computational requirements in multicanonical 
simulations dramatically, when compared with the full range calculations in Sec. (|4.3| ); 
using random walk arguments, the reduction factor is 0.15^ ~ 0.02; in practice, the 
reduction is even stronger than this. 

In Fig. ^we show the actual order parameter distributions P(0^y), measured from the 
aX2 = 3, 803 lattice, and reweig hted to supercooling values 6m^/X2 = 0.00077-0.00122. 
The larger the supercooling, the smaller the critical droplet, and less suppressed the 
droplet probability is. We always use lattice units for 0^^; since this is an "internal" 



quantity in the expression for the rate ( p. 17] ), there is no need to convert it to continuum 
units. 

We calculate the critical droplet free energy through Fc ~ — ln(P(0^^ (^)/P(0^y J), 
where P(0^y J is the height of the symmetric phase peak of P(0^y). Note that this 
ratio is not P^, which appears in the rate equation ( p^.l7D ; we chose to use this ratio 
because it is dimensionless, whereas the probability factor in (|2.17|) has dimensions of 



[^av]~^- Fig. we show Fc using A2C1 = 3 lattices of volumes 603-803. For large 
values of supercooling 6m'^, where the critical droplet is small, we do not observe any 
significant finite volume effects. However, when (5m^ is smaller, the critical droplet is 
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Figure 6: The probability distribution -P(0av) various values of supercooling (5m^, 
measured from aX^ = 3, 80^ lattice. Here is in lattice units, Eq. (|3.7|) . The critical 
droplet value q at each Sm"^ is defined to be the location of the local minimum of 

large and, in small lattices, the droplet feels the proximity of its periodic copies. For 
our smallest (60^) lattice, this decreases the droplet free energy significantly already 
at modest droplet sizes. The curves shown in Fig. ^ correspond to droplets which are 
well within the maximum size determined by the 15% rule discussed in Sec. (|2.1| ). This 
finite size effect is caused by the finite thickness of the droplet wall, and it should 
vanish exponentially as the lattice size is increased. Indeed, at large lattices we need 
to go to very large droplets in order to observe any significant finite size effects. 

On the other hand, the lattice spacing effects are substantial for the values of aX2 
used: in Fig. | we compare Fc from lattices of similar physical size, but aX2 = 2, 2.5 
and 3. The larger lattice spacings aX2 = 3.5 and 4 have much larger finite a effects and 
we discard these in our analysis. For a fixed value of supercooling Sm^, Fc decreases as 
the lattice spacing becomes smaller. The leading order errors are expected to be 0{a^)] 
indeed, a good fit to the 3 curves is obtained with a 2 parameter fit Cq + C3(aA2)^, done 
independently for each 6m'^. The resulting continuum curve is shown in Fig. |[ The 
statistical errors in the Fc curves are ~ ±0.5-1.4, and in the continuum extrapolation 
^ ±2; this is shown in the figure. 
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Figure 7: The critical droplet free energy Fc-, as a function of supercooling 5m?, 
measured from Xoa = 3 lattices of size 60^-80^. 
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Figure 8: The critical droplet free energies measured from = 3, volume 80^; 

= 2.5, volume 100'^, and = 2, volume 120^ lattices. The continuum (a = 0) 
curve is shown with an error band. 
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5.2 Real time evolution (i): {\A(pl^/At\) 

The factor {\A(j)l^/ At\) is to be evaluated at the critical droplet value = 0avC- 
is easy enough to measure from numerical simulations, but it can be also calculated 
analytically: from the equations of motion ( p.8|) we see that 

At m 



Y.^a<iya + o{{Aty). (5.1) 



x,a 



Here all the fields and At = (At)cont./'2 are in dimensionless lattice units. In thermal 
equilibrium, the momenta na{x) have Gaussian random variable distribution (of width 
one in our normalization). This implies that {7iaix)7ib{y)) = Sa,bSx,y, and any expecta- 
tion values involving products of vr's and 0's factorize: {f{(j))g{7!-)) = {f{(f))){g{7i)). We 
need to determine (|A0^y/At|) at the special value = 0avC- Using this constraint, 

Ex,a0a(^) = ^^0av,C> ^C obtalu 



2 N 2\ A /l^2 



For any fixed x the product 7ra{x)(f)a{x) has a non-Gaussian distribution, but the sum 
which appears in Eq. ( p.l[ ) is Gaussian, due to the global constraint. Because Gaussian 
distributions satisfy (x^) = 7r(|x|)^/2, we obtain the result 



A0 



At 




(5.3) 



where we have converted back to continuum time, to remind us of the correct scaling 
in the continuum limit. Since this result depends only on the average distribution of 
momenta vr, it is independent of the magnitude of the noise in Eqs. (|3.8|) . 

5.3 Real time evolution (ii): (d) 

The final contribution to the rate comes from (d) = (5tunnei/A^crossings)- This requires 
the evaluation of full real time trajectories through an ensemble of critical droplets. 
Following the procedure outlined in Sec. (|2.3|), we do this as follows: 



(1) First, we choose an initial configuration 0q(x) from a thermal distribution, but with 
the order parameter restricted to a narrow interval around the critical droplet value: 
\4>'i-y — (p'iy c\ ^ ^/2- These are straightforward to generate with either a standard Monte 
Carlo simulation with the restriction built into the update, or by choosing them from 
a multicanonical run. 
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Figure 9: Two trajectories measured from = 3, 80^ lattice, at supercooling 
Sm'^/Xl = 0.00094 and with noise magnitude = 27 At = 0.00125. The trajectories 
are evaluated forward and backwards in time, starting from configurations at t = 0. 
Trajectory A crosses 0avC times, B 5 times; the contribution to d from A is 0/12, 
and from B 1/5. The dashed horizontal line is the symmetric phase expectation value; 
and (f)l^ is shown in lattice units. 

(2) We assign initial momenta vrj](x) to this configuration, drawn from a thermal dis- 
tribution. In our case this means that we choose iTaix) from a Gaussian distribution 
with width a^(7r^) = 1 (or (tt^) = 1 in lattice units). 

(3) The configuration is then evolved in time until the order parameter (pl^ reaches 
'^av,s ^1^,3^ symmetric or broken phase "cut-off" values. We always use the 
timestep At = 0.05a in our real time runs. We then return to the initial configuration 
(j)^, invert the initial momenta 7r° —>■ —tt^, and evolve the system again until we reach 
4>'iv,s 0av,B- This latter run is interpreted as a run backwards in time; by gluing the 
backwards and forward half-trajectories together at the starting point, we obtain a full 
trajectory (pl^^S or B) — > (f)l-^{S or B). If the trajectory tunnels, i.e. if the backward 
and forward evolution ends are on different sides of q, it contributes to the tunneling 
rate (5tunnei = !)• After counting the number of times the trajectory crosses (f>1^cy ^® 
obtain its contribution to (d). 

An example of the trajectories is shown in Fig. |^. About 40% of the trajectories we 
evaluate are tunneling trajectories. This means that the order parameter resolves the 
critical droplet quite well, and configurations at = q indeed are a good sample 
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Figure 10: The dynamical factor d plotted against 5rn? , measured from aX^ = 3, 80^ 
lattice. 

of critical droplet configurations (see the discussion in Sec. ^^) . We measure the order 
parameter after every timestep. On the scale of the plot the evolution of (pl^ looks 
rough, but at the level of individual timesteps it is quite smooth. This is due to the 
relatively small amplitude of the noise in the equations of motion. 

In Table |I] we show the lattices and Srri^ values where we evaluate the trajectories. 
In Fig. |iy we study how d depends on the degree of supercooling 5m^. Remember 
that small values of correspond to large critical droplets. At each value of Sm"^ 
we evaluated 35-83 full trajectories. The value of d appears to be remarkably stable 
throughout the range of studied; the free energy of the critical droplet varies by 
~ gioo Q^gj, ^]-^jg Ya,nge. Naturally, d still depends on 5m^; for example, when —>■ 
one expects d — 0, simply because the larger the droplet is, the slower it evolves. 
However, the dependence on Sm^ is not visible within our statistical errors. The large 
values of the errors are caused by the large variation in the number of crossings between 
tunneling trajectories: some trajectories have only a few crossings, whereas some have 
more than 50. The errors of d are still sufficiently small for an accurate calculation of 
the the nucleation rate — since the rate is of order e~^°°, a factor of 2 means little in 
the final result. 

The distribution of the trajectories as a function of the number of crossings appears 
to decrease roughly exponentially as the number of crossings increases. However, when 
the number of crossings is small, the trajectories with an even number of crossings — 
which do not lead to tunneling — occur more frequently than those with an odd number 
of crossings. This is probably a consequence of the order parameter fiuctuations in the 
bulk phases, discussed in Sec. 
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The lattice spacing dependence of d is small: in addition to the aX2 = 3 results 
in Fig. we have measured it using lattice spacings aX2 = 2.5 and a\2 = 4. At 
aX2 = 2.5, with a lattice of size 100^ and supercooling 5m? /\\ = 0.00094, the result is 
(d) = 0.049 ± 0.022; and at aAs = 4, volume 60^ and supercooling 6m'^/Xl = 0.001156 
the result is (d) = 0.064 ± 0.018. Thus, no significant lattice spacing effect is seen. 
Note however that on extremely fine lattices, (|A0^^/At|) oc ~ a~^^^ diverges 

when expressed in physical units (if we keep At smaller than the lattice spacing), and 
d must correspondingly go to zero as a^/^. 

5.4 Effect of the noise 

Let us study how the rate is affected by the magnitude of the noise + damping terms, 
parametrized by 7 in the equations of motion ( p.8|) . Since different levels of noise 
correspond to different dynamical evolution (also in the continuum), there is no reason 
why the results could not have significant dependence on it. As mentioned before, 
(|A0^y/At|) is independent of the noise, as is P^cpl^), so it is sufficient to study how d 
depends on 7. 

What kind of behavior can we expect? Naively, neglecting the interactions and the 
mass (Ai = A2 = m'^ = 0) and the noise term ^ in the equations of motion, the evolution 
of a mode of wave vector k obeys — i'yuj — = 0. If we have small 7 < 2k, we see 
that |ti;| = k independent of the value of 7. On the other hand, if 7 > 2k, the evolution 
becomes overdamped, and if 7 is very large u ~ ik'^/j. Boldly extrapolating these 
simple arguments to the physical evolution of the critical droplet, we can expect that 
at small 7 the rate is approximately constant as 7 increases, whereas above a threshold 
value it should behave as I/7. 

In Fig. |TT]we show the measured value of d against e = (1— e"^"^^*)^/^. Except for the 
points near e = 0, we indeed observe roughly the expected behaviour, and, surprisingly, 
even the threshold scale 7 ~ 2A; 27r/(droplet size) ~ 0.25 (or e ~ 0.16 in Fig. |lT]) 
is reproduced by the data. At e = the evolution is Hamiltonian, and, as discussed 
in Sec. (|3.3| ), the finite volume causes additional complications: the growing/shrinking 
droplet releases/absorbs latent heat, which rapidly equilibrates throughout the system. 
Since the Hamiltonian evolution conserves total energy, on a finite volume this causes 
slight heating/cooling of the system. This, in turn, tends to stabilize the droplet, 
strongly reducing the tunneling rate. Indeed, in a small enough volume the critical 
droplet may not decay at all. 

The addition of the noise to the equations of motion thermalizes the system effec- 
tively. A natural amplitude for noise is 7 = l/L, which thermalizes the system in 
the same timescale as waves propagate through it. This corresponds to e = 0.035 in 
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Figure 11: The dependence of d on the noise magnitude e = (1 — e"^"^^*)^/^, measured 
on an a\2 = 3, 80^ lattice at supercoohng Sm'^/Xl = 0.00094. The point at e = 
corresponds to Hamiltonian evolution, and the "heat bath" result to completely noisy 
evolution e = 1, shifted, for clarity, to smaller e. 

Fig. |TT|. Indeed, from this point up to £ ~ 0.16 d does not vary much, and certainly 
the variation is not significant for the final tunneling rate calculation. Thus, we can 
assume that in a very large volume, the Hamiltonian evolution would also give a value 
for d at this level or slightly larger. 

5.5 Nucleation rate 

Finally, let us pull together the ingredients discussed above and calculate the rate F/A^ 
with Eq. ( p.lTp . First, note that since the dimension of Pq is [0aviat]~^ = ['^^L cont.]""^) 
it does not have a good continuum limit by itself. In order to cancel the factor, it is 
convenient to multiply Pc/iVXl) by (aA2) before extrapolation, and correspondingly 
divide (|A02^/At|)(d) by it. (This is not an issue for -P(0av c)/-^('i^av sym) used in 
Sec. pTT| , since it is dimensionless.) 

In Fig. |1^ we show the nucleation rate F/Ag from aX2 = 3 lattices using 7 = 0.0125, 
where we have the most extensive set of data. In the "probability" curve we have set 
the product of the dynamical factors l/(2aA2)(|A0^y/At|)(d) equal to A2 (this gives 
correct dimensions), and only the probability contributes to the rate. In the "full 
rate" curve we include the correct value of {\A(j)1^/ At\) from Eq. (|5.3| ), and d from 
Fig. |10|, by substituting it with its average and making a conservative error estimate, 
d = 0.08 ± 0.04. In the final error estimate we also take into account the errors of PA. 
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Figure 12: The rate of droplet nucleation from a\2 = 3 lattices. The full rate 

curve contains all of the terms in Eq. ( |2.17D , whereas the "probability" curve contains 
only the contribution Pc/{V^2)- The thin wall curve is calculated using the surface 
tension and the latent heat measured from aX2 = 3 simulations. 

The inclusion of the dynamical factors reduce the rate by a factor of ~ 10^. However, 
we emphasize that only the full rate is independent on the choice of observables. For 
example, if we simply switch from an intensive to an extensive order parameter, Pq, 
which is proportional to l/^av would decrease by a factor oc l/V; this is compensated 
by a corresponding increase in (|A02^/At|), so that the full rate remains invariant. 

Next, let us consider the continuum limit extrapolation of the rate F/A^. The ex- 
trapolation of the probabilistic factor PQ{a\2)/(y\l) alone proceeds as in Fig. ^; we 
fit an ansatz of form cq + C2a^ independently at each (5m^/ to the lattice results from 
aX2 = 2, 2.5 and 3. The resulting curve is shown in Fig. |13[ The remaining contri- 
bution to the rate is given by the dynamical factor (aA2)""'^(| A0^.^/At|)(d), which we 
assume to be constant as the lattice spacing is changed but the physical volume is kept 
constant. This quantity must have a continuum limit unless the dynamics posess some 
UV pathology. The simulation results from different lattice spacings do not show any 
lattice spacing dependence within statistical errors; however, errors are too large for a 
real continuum extrapolation. 

The result of the extrapolation is shown in Fig. |T3|. The errors of the log of the rate 
are ~ ±4, including the errors coming from the extrapolation of P^ and the estimated 
error from the extrapolation of d. The error in Pq dominates the uncertainty. 
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Figure 13: The nucleation rate as in Fig. but extrapolated to the continuum limit. 
For comparison, we also show the rate calculated using the thin wall approximation 
and leading order and next-to-leading order perturbation theory. The curves labelled 
K = 0.5 and 1.5 show the scale sensitivity of the NLO calculation, see Sec. [Oj 



5.6 Droplet cross-sections 

It is illuminating to take a closer look at the structure of the droplet configurations. 
Since the droplets are large objects in terms of correlation lengths, it is possible to 
perform suitable averaging over lattice configurations in order to measure the size and 
shape of critical bubbles at fixed supercoohng. We do the averaging as follows: 

• Generate a number of configurations with the constrained order parameter — 
^avcl ^ ^1 where e is a small number. These configurations contain a droplet of 
a fixed volume, as discussed in Sec. (^]l|). 

• Determine the center of the droplet. We do this by averaging the order pa- 
rameter over z coordinates, and taking the lowest non-trivial Fourer mode 
^ = • The value xq = L/n arctan(Im74/Re74) gives now the 
point around which the configuration is maximally symmetric. This is repeated 
for y and z. 

• Shift the origin to the center of the droplet, and determine the average order 
parameter as a function of the radius, 0^(r). Averaging over the configurations 
we obtain the average droplet cross-section (0^(r)). 
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Figure 14: The droplet cross-sections at 3 different supecoolings from aX2 = 3, volume 
80^ lattices. The expectation value is shown in lattice units. The dotted line shows 
the broken phase at the critical value of m^, the solid lines at the supercooling 
(5m^. Only at the smallest supercooling does the core of the droplet approach the 
broken phase expectation value. 

We determined the droplet cross-sections for several different at aX2 = 3 lattices 
of size 80^, and in Fig. we show the results for supercooling 5m? /\\ = 0.00124, 
0.00094 and 0.00077. From Fig. |^ we see that these correspond to droplet free energies 
from 75 to 180. The droplet wall is very thick, and especially for smaller droplets 
there is no broken phase core to speak of.[^ Thus, the thin wall approximation for 
the nucleation rate cannot be expected to work well. The expectation value (0av) 
larger at the center of the droplet than the broken phase expectation value at m^; 
this is naturally because the broken phase expectation value increases rapidly as is 
decreased. Only for the largest droplet (largest m^) does the core approach the broken 
phase expectation value. 

One can note that although the bulk of the droplets fits well on the lattice, at a 
smaller volume the exponentially decreasing tail of the fields outside the droplet wall 
would feel the lattice size. This decreases the free energy of the droplets, and it is 
no surprise that the 60^ lattice gives smaller droplet free energy at large droplets, see 

^^One should be a little cautious of interpreting the interface width shown in Fig. however, 
since an interface in 3 dimensions generically exhibits logarithmically large roughening. Locally the 
interface may be thinner. 
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Fig. 0. 

6 Comparison with other nucleation rate calculations 
6.1 Thin wall approximation 

The simplest semi-analytical estimate for the nucleation rate can be obtained using the 
thin wall approximation for the critical droplet free energy Eq. (|2.5|): 



where we use lattice determinations of the surface tension a and latent heat ^. Thus, 
this method is not fully analytical; however, normally the determination of these values 
is an order of magnitude easier task than the full nucleation rate computation, and the 
thin wall approximation has often been used to include non-perturbative input in the 
nucleation rate calculations (see, for example, p3| ). 



The free energy is converted to a rate by multiplying it with suitable mass scale 

r = m"^ exp — Ftvk ■ (6.2) 

We use here the mass m = O.O5A2, which is close to the symmetric phase perturbative 
mass; the precise value is not significant. The thin wall approximation assumes that 
the thickness of the droplet wall is negligible and that there are no contributions to 
the free energy due to the curvature of the walls. It becomes valid in the limit where 
the droplet radius is much larger than the thickness of the droplet wall; this condition 
is not usually well met in practice (see Sec. ^]6|). 



In Fig. |T2| we show the thin wall result using the values of a and i obtained from 
aX2 = 3 simulations. In the range of interest the thin wall rate is higher by ~ e^^. 
Large as the difference is, the thin-wall approximation is not as bad as one may first 
think: in many physical processes the system is cooled down at a very slow (constant) 
rate, and thus, the nucleation rate is changing with time: dT/dt > 0. Since F depends 
exponentially on 5m^, the transition occurs very rapidly when the nucleation rate 
reaches some critical value, typically F^ x l^[t — tc) ~ 1, where I is some relevant length 
scale. 

Let us take F/A2 = e~^°° as our reference value. The lattice results from A2 = 3 give 
a supercooling value 6m'^/\l = 0.00132(2), whereas the thin wall result is 0.00096(4), 
which is "only" ~ 30% smaller than the correct value. The errors of the thin wall 
results come from the errors of the determination of a and £. 
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In Fig. |T3 we compare the lattice results to the thin wall ones at the continuum limit. 
Using again the reference value of Fc/Al = e~^°°, the lattice result for the supercooling 
is 5m'^/\l = 0.00094(3), and for the thin wall calculation 0.00063(7), again a 30% 
difference. Thus, the thin wall calculation gives a rather good order of magnitude 
estimate for the nucleation rate. This was also seen in the SU(2)-Higgs calculation in 



Ref. 



6.2 Perturbation theory 

In the broken phase of the cubic anisotropy model (and with our choice of A1/A2) the 
field component which acquires a non-zero expectation value (0i, say) is much lighter 
than the other field component (02)- Due to this mass hierarchy, it turns out that 
we can calculate an effective potential (or action) for the 0i field by perturbatively 
integrating over 02- This effective potential can be used to calculate the leading con- 
tributions to the critical droplet free energy by finding the classical "bounce" solution 
to the equations of motion |3^, |35|, More concretely, assuming that the droplet 
is spherically symmetric and centered at r = 0, we want to find a configuration v{r) 
which is a saddle point of the action 



S{m^) = An drr^ 
Jo 



(6.3) 



Here v is the expectation value of the scalar field. The boundary conditions are f (00) = 
0, drv{0) = 0. 

The dimensionless expansion parameter is proportional to the ratio of the two cou- 
plings A1/A2. This implies that the convergence becomes better the stronger the first 
order transition is (see Sec. For our parameters the expansion parameter is 1/8 <^ 1, 
offering formally a good convergence parameter. 

Without loss of generality we can choose the transition to happen in the direction of 
the 01 field, and we shift the fields (0i, 02) {v + 02). The tree (mean field) level 
potential is 

Vo{v) = ^mh' + ^X,v^ (6.4) 



and the 1-loop potential is 



where 



Vi{v) = - — {ml + ml), (6.5) 

iZTT 



1 1 
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In the broken phase, for m? close to the transition value, inspection of the potential 
gives the following power counting rules: 

v"^ ~ Xl/Xj and ~ A^/Ai . (6.7) 

The leading contribution from the 1-loop term is of the same order as the tree 
level terms; it is just this third order term which makes the transition first order. On 
the other hand, the contribution from mf is suppressed by a factor oc (Ai/A2)^''^ in 
comparison to the leading order. This is actually subleading when compared with 
the leading 2-loop contributions; furthermore, at this order one also gets contributions 
from the resummation of the (pi propagator, which render mi formally imaginary (easily 
seen by considering Vi{v), and equating this to ml). It actually turns out that up to 
next-to-leading order in A1/A2 we do not have to consider the fluctuations of (pi at all, 
justifying the use of the classical "bounce" solution at this order. 

At the leading order we can also neglect in the expression for m"^, and write the 
potential as 

Vi^oiv) = Vo{v) - ^^xTv\ (6.8) 

This potential gives a very strong first order transition, much stronger than the lattice 
calculations indicate, as can be seen by comparing the broken phase values v"^ and the 
surface tension a in Table 0. We calculate the perturbative surface tension from the 
integral 



a 



^broken 



dv^2V{v) . (6.9) 

' ''^symm. 

The droplet free energy (at some value of supercooling) can be solved from Eq. ( |6.3|) 
numerically using the well-known undershoot /overshoot method [^. (It can be solved 
analytically in some limiting cases, see, for example, |3^.) As in the thin wall droplet 
case, we convert the droplet free energy to the nucleation rate using T = m^exp(— F), 
where we use the mass scale m = O.O5A2 (the result is insensitive to the precise value 
of m). Solving for the value of supercooling which gives the reference rate r/A2 = 
exp— 100, we obtain the supercooling Sm"^ = 0.00234A2, which is more than 100% 



larger than the result from the lattice simulations. Fig. |T3| and Table |[ 

The next-to-leading order corrections to the potential are suppressed relative to the 
leading terms by A1/A2. They arise from in Eq. ( |6.5|) , and from the leading two-loop 
contribution: 



1 , v^X^ 



V,{v) = Vo{v) - —mi - 



3 ^'2 



1 /i^ 1 



(6.10) 



127r ^ 4(47r)2 

The last term above is the two-loop contribution; it arises from the logarithmically di- 
vergent "sunset" diagram with two (p2 -propagators and one 0i-propagator, see Fig. |15|. 
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Figure 15: Left: 2-loop "sunset" diagram used in the effective potential calculation. 
Right: the diagram used in the calculation of the wave function correction. The dotted 
lines are the "heavy" 02 -propagators, and the solid lines correspond to the "light" 01 
propagators. The coupling constant at the vertices is X2V. 
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Table 2: Comparison of the perturbative, lattice and thin wall results. The supercooling 
has been evaluated at the point where the rate T /\\ = exp(— 100). 



At this order we can neglect the masses from the propagators in the 2-loop diagram, 
and also the other 2-loop contributions (sunsets and figure-8 diagrams), which are 
suppressed by {Xi/\2Y^^ or more. 

In Eq. ( 6.101) /i is the scale factor arising from dimensional regularization in d = 
3 — 2e dimensions. As usual, a logarithmic divergence has been absorbed in a mass 
counterterm, which defines the scale dependent running mass 

We substitute with m|, everywhere where it appears in the potential. 

The dependence of V{v) on the scale fi is formally of higher order than A1/A2; thus, 
we are free to choose fi so that the effect of the large logarithms in the potential are 
minimized. A common prescription is to fix fi so that the logarithm vanishes in the 
broken phase. Naturally, this does not suppress the logarithms in the symmetric phase 
at all. A better overall cancellation of the logs is achieved by using a f-dependent scale 



/i = 2k J X2v'^/2, where k is a constant of order unity. However, this simple prescription 



is not complete, as discussed in Sec. 5 of Ref. [37|. Simplifying the renormalization 



arguments in the above reference a bit, this follows from the fact that the value of the 
effective potential is not a physical quantity, only the differences V{vi) — V{v2) are. 
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and the absolute value of the potential can have unphysical /i-dependence. (Indeed, in 
dimensional regularization the value of the potential also at v = depends on the scale 
/i.) Thus, if we change the renormalization point /x as we change v, the "normalization" 
of the potential changes, and it is not possible to directly compare V{vi, fi{vi)) and 
V{v2, fJ,{v2)) if Vi and V2 are far away from each other. 

The most straightforward solution to this problem is to consider the slope of the 
potential dV{v)/dv instead of the potential itself. With this quantity we can write 
a renormalization group improved expression for the next-to-leading order effective 
potential p7[] : 

dV2{v,fi) 



dv 



/i 



(6.12) 



The derivative in the above expression is to be evaluated at fixed /i, and only after 
taking the derivative /i is set to its optimized value. 

In Table |^ we show f broken ^^'^ calculated from potential (|6.12|) , using values 
K = 0.5 and 1.5. The transition is weaker than the results from the lattice simulations; 
the broken phase expectation value of ^broken is 10-20% smaller, but the surface tension 
is around 40% smaller. Likewise, the supercooling 5w? at rate V = A2exp(— 100) is 
around 50% smaller than the lattice results. This discrepancy is of similar magnitude to 
the difference between perturbative and 2-loop results in the nucleation rate in SU(2)- 
Higgs theory [|1^]; however, in contrast to the cubic anisotropy model, the supercooling 
at a fixed value of the rate was seen to be larger in the perturbative analysis than on 
the lattice. 

The scale dependence of the NLO nucleation rate is quite large in Table ^ the value 
of the supercooling increases by 75% when k is decreased from 1.5 to 0.5. This suggests 
that higher orders would affect the results significantly — this is, of course, already 
indicated by the fact that the NLO results do not agree well with the leading order 
results, or with the lattice. The large scale dependence is not a feature of the procedure 
to set fi used in Eq. ( |6.12| ): we can, for example, use the 2-loop potential ( |6.10| ) and 
just fix the scale /i so that the 2-loop contribution vanishes in the broken phase. If we 
vary by a factor of ~ 50%, the strength of the transition changes about as much as 
when we vary k in Eq. ( p.l2| ). Moreover, this method tends to give a somewhat weaker 
transition than Eq. ( 6.12|) . 

In principle, finding the extremal critical droplet solution of the bounce action 
Eq. (|67^ ) with the NLO potential ( |6.12D does not give all of the 0{Xi/X2) contri- 
butions. At this order one should take into account also the wave function correction 
to the kinetic energy in the bounce action by modifying the second derivative term: 

{Vvf ^ (1 + Z^){\/v)\ (6.13) 
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The leading contribution to can be calculated by evaluating the diagram on the 
right in Fig. |15| and expanding the result to second order in the external momentum 
(see, for example, Ref. This calculation is justified because v (or 0i) varies on 



length scales ~ t^i^-, which is much longer than the inverse mass of 02, ^^2^- The 
result of the calculation is 

which is indeed oc A1/A2 according to the power counting rules, Eq. (|6.71 ). However, the 
inclusion of the wave function correction modifies the results by less than 1%, which is 
a negligible effect considering the accuracy of the NLO calculation. This is in contrast 
to the SU(2) gauge + Higgs theory, where the wave function normalization was seen to 
improve the results dramatically. The reason for the smaller effect in this case is the 
very small numerical multiplicative factor in Eq. ( |6.14] ). 

We can conclude that the NLO perturbative analysis at A1/A2 = 1/8 describes the 
phase transition qualitatively correctly, but the physical quantities can be a factor 
of 2-3 off the correct ones. Naturally, for A1/A2 closer to unity the NLO analysis 
would be even less accurate. The overall accuracy is comparable to the perturbative 
treatment of the SU(2) gauge + Higgs theory [|19[. Trying to improve the perturbative 
treatment would require the calculation of contributions proportional to (Ai/A2)^^^. 
At this order the light 0i loops start to contribute, and, as mentioned above, 0i 
propagators require resummation, which makes mi (and the 1-loop contribution in 
( |6.5| )) formally imaginary in the region where the second derivative of the effective 
potential is negative. Thus, it is not consistent to solve the effective potential in a 
constant background field as above; the correct treatment would require the calculation 
of the fiuctuation determinant around a classical droplet solution. The simultaneous 
consistent (as a power series of A1/A2, say) evaluation of the fiuctuation determinant 
and the effective potential is a very non-trivial problem, and we do not attempt it 
here.|^ 

6.3 Comparison with the coarse-grained effective action approach 

The nucleation rate in the cubic anisotropy model was studied by Strumia and Tetradis 



using a coarse-graining approach [|T4|. First one calculates a coarse-grained effective 
action, valid for momentum scales smaller than a chosen scale k, where k<m, some 
physically relevant mass scale (for example, the mass of the light 0i field in the broken 



^^The fiuctuation determinant iias been recently performed by several authors [^9[ 14, ^ in 
various physical systems. However, these studies did not consider the simultaneous order-by-order 
evaluation of the perturbative effective potential and the fluctuation determinant. 
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phase, if the transition happens along the direction of (pi). Formally, the effective action 
is found by integrating out momenta p > k. In the second step, the nucleation rate is 
calculated using the effective action and Langer's method: one finds the classical bounce 
solution (as in Sec. 6^) with action Sk, around which the fluctuation determinant 



can be evaluated. The scale k acts automatically as the ultraviolet regulator for the 
fluctuation determinant. The full rate is then estimated to be 

r = AkC-^'^ . (6.15) 

Since the physics does not care about the artificial scale k, the final result should be 
independent of k if both halves of the calculation are under control. 

The authors of Ref. |n] find the effective action using the renormalization group 



flow method of Wetterich they start from the bare action, Eq. ( |3.2| ), defined at 



some ultraviolet scale k^ ^ fc, and use the flow equations to evolve the action down to 
scale k. During the evolution the action develops so that the resulting effective action 
has a first order transition already at "tree level" . It still contains both of the original 
fields. The evolution generates very complicated (non-local) effective actions, and in 
practical calculations the actions have to be truncated to a simple ansatz. For details, 
see Ref. and references therein. 



In Ref. the nucleation rate is calculated at coupling constant ratios A1/A2 = 0.3 
and 0.15 (in their notation, Xko/gko = 0.1 and 0.05, defined at the ultraviolet scale ko.). 
The latter is close enough of the value we use here A1/A2 = 1/8 to make qualitative 
comparison possible. The results show a very strong dependence on the coarse graining 
scale k: for a fixed supercooling, the nucleation rate varies from InF/m^ ~ —250 to 
— 130, when k/mi is changed from 0.6 to 0.9. Here mi is the mass of the light 0i 
field in the broken phase. This implies that accuracy is lost during the calculation. 
Furthermore, the leading contribution to the rate does not come from exp{—Sk), but 
from the fluctuation determinant Aj.. This is partly caused by the fact that the mass 
of the 02 field in the broken phase is much heavier than the coarse graining scale k, 
and gives a large contribution to the determinant. Indeed, one can argue that a coarse 
grained effective action where the heavy field is completely integrated out, as was done 
in Sec. |6.2| , may offer a better starting point for the Langer method.]^ 

We should emphasize that the approach just described makes two separate approx- 
imations. The first is the truncation of the infrared effective action, made in the the 



^''This procedure has been suggested by Weinberg Q: one obtains an effective action for the hght 
fields by integrating over heavy fields perturbatively. Only the light field fluctuation determinant is 



evaluated. However, as discussed in Sec. 6.2, this is very difficult to implement as a correct order-by- 
order perturbative expansion. 
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renormalization group flow part of tlie calculation as the coarse graining scale is run 
down to k. This approximation means that some ultraviolet physics is lost or incor- 
rectly incorporated in the coarse grained effective theory. The second approximation 
is in the calculation of the nucleation rate within the coarse grained theory itself; com- 
puting Gaussian fluctuations about a classical saddle point corresponds to carrying 
perturbation theory to one loop, neglecting interactions between infrared fluctuations. 
This is clearly not warranted if the infrared behavior is strongly coupled, for instance. 

Since the dominant contribution does not come from the saddle point but from 
the fluctuation determinant, hanger's method as a saddle point expansion fails in 
this calculation. However, our lattice simulations clearly indicate that a well-defined 
saddle point exists in the phase space. While the simulations make no distinction 
between the classical solution and the fiuctuation determinant (everything is contained 
in the droplet free energy), the configurations at the saddle point consist of well- 
defined critical droplets. Thus, the physical picture given by hanger's theory is valid; 
the problem in (semi) analytical calculations is to find the correct effective description. 



7 Conclusion 

In this paper we have used a novel technique for studying the critical droplet nucleation 
rate in first order phase transitions using non-perturbative lattice simulations. The 
method is readily applicable to exponentially small nucleation rates, which are often of 
relevance in physically interesting transitions. The technique consists of two separate 
stages: (1) the Monte Carlo evaluation of the nucleation barrier, using multicanonical 
methods, and (2) the correct treatment of the microscopic dynamics of the nucleation 
process with real time simulations. The first stage can be considered as a generalization 
of hanger's theory of nucleation : we have substituted the approximate saddle point 
expansion with a nonperturbative Monte Carlo calculation. The second step goes 
partly beyond hanger's formalism. This method is readily applicable to any theory 
which is amenable to the lattice treatment in the first place, and, within the context of 
the lattice, the technique is exact up to exponentially small and in practice negligible 
corrections. This method has also been applied to the nucleation rate in SU(2)-Higgs 
theory Due to the multicanonical methods used, the statistical errors in the final 



answer tend to be automatically small — usually much smaller than any uncertainty 
in analytical approaches. 

In this work we have applied the method to the phase transition in the cubic 
anisotropy model, which is a formally simple field theory with two scalar fields and 
a radiatively induced first order phase transition, the strength of which can be ad- 
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justed continuously. Because of the formal simplicity the theory is well suited for both 
analytical and numerical analyses. Since the transition is radiatively induced, Langer's 
theory has to be applied with great care in analytical calculations. 

We compare the nucleation rate obtained from the lattice simulation to customary 
analytical or semianalytical approaches. We find that the straightforward application 
of perturbation theory up to next-to-leading order in the relevant expansion parameter 
is not accurate: for a fixed value of the nucleation rate, the corresponding supercooling 
is roughly a factor 2 off the correct value, even though the scalar field expectation 
value is correct to within 15%. To go beyond the next-to-leading order would require 
the evaluation of the fluctuation determinant, and we did not attempt it in this work. 
On the other hand, a thin wall estimate, using non-perturbative input, is off only 
by 30%. The nucleation rate for a fixed supercooling is off even more than these 
numbers. This behaviour is strikingly similar to that observed in SU(2)-Higgs theory. 
The nucleation rate in the cubic anisotropy model has also been studied in Ref. [1T4| , 



using an approximate coarse grained effective potential as a starting point for Langer's 
procedure. The results were seen to depend very strongly on the coarse graining scale, 
making quantitative comparison impossible. 

The results show that, for a rough-and-ready estimate of the nucleation rate, the 
thin wall approximation is acceptable, provided that one uses non-perturbatively deter- 
mined surface tension and latent heat as an input. These are much easier to determine 
on the lattice than the full nucleation rate. On the other hand, the purely perturbative 
analysis shows weak convergence. If high accuracy is required, one has to resort to 
numerical evaluation. 
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A Renormalization of the lattice theory 

This appendix details the matching calculation which eliminates O(a^) errors in the 
lattice theory. This is done by computing a set of correlation functions in the lattice 
and continuum theories, at one and two loops, and adjusting the lattice couplings so 
the results coincide. The required correlation functions are the two and four point 
functions at zero momentum, the leading momentum dependence of the two point 
function, an insertion of (0^) on a zero momentum line, and the vacuum value of the 
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{(fp') operator. It is necessary only to find the difference between lattice and continuum 
values, to perform the matching; this difference is IR finite for all graphs we need. 
Each loop order eliminates errors at one higher power of a because graphs become 
more UV convergent by one power per loop order in this 3 dimensional theory. For a 
more thorough discussion of how the matching calculation works, see . 



A.l results 

A one loop lattice-continuum matching (renormalization) calculation will determine 
the 0(a) contributions to the quantities Z^, 8\ and Z^, and will find the 0{l/a) con- 
tributions to 5m? and The required graphs are shown in Fig. |16|. Evaluating the 
graphs requires choosing a lattice Laplacian. We consider two choices; the unimproved 
Laplacian 

VU{x) = -60(x) + 5:(0(x + 0(x - i)) , (1.1) 

i 

and the improved Laplacian, which we actually use: 
ViV(x) = -y 0(a;) + ^ E {^i^ + + - - ^ E {^i^ + 20 + (i>{x - 2i)) . (1.2) 

i i 

We write the Fourier transforms of these choices as 

= E(2-2cos(A;i)), (1.3) 

i 

= E - ^ ^°^(^*) + I cos(2A;,)) . (1.4) 



The evaluation of the one loop graphs requires two integrals: 



Att Jbz (27r)3 A;2 
^ _ f d^k 1 [ d^k 1 



47r Jbz (27r)3 {k^y J^s {2-nf 



:i.5) 
:i.6) 



Here we use the shorthand BZ to mean k lies within the Brillouin zone, meaning each 
ki G [— TT, tt]. The integrals which determine ^ are each IR singular and some regulation 
is implied, for instance adding m? to both k"^ and k"^ and taking the limit as m? 0. 
The numerical values of the integrals are 

Su = 3.17591153562522, Sj = 2.75238391130752 , 

= 0.152859324966101, = -0-083647053040968 . (1.7) 
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(c) (d) 

Figure 16: The one loop graphs needed for the renormahzation. A cross represents a 
(f^ insertion. 

Note that the sign of ^ depends on whether we use an improved lattice Laplacian. This 
is possible because ^ represents the difference of a graph between lattice and continuum 
theories. The lattice contribution is larger inside the Brillouin zone, but the continuum 
integral receives contributions from outside the zone as well; the sign depends on which 
effect is larger. 

At one loop the renormalizations are 

«M, - (1.8) 

«2,u = (A,A2 + 2Al) , (1.9) 
Z<,,u-1 = 0, (1.10) 
Z„u,-1 = Q^. + ^A,)!^. (1.11) 

Sml = -(iA, + lA.)^. (1.12) 

= 2^-2m^4. (1.13) 

Note that, if Ai = A2, then 5\i = 6X2; and similarly if Ai = 3A2, then 6\i = 36X2. 
Therefore the decoupling and 0(2) symmetric versions of the theory are preserved 
under renormahzation. 

It makes no sense to carry the matching to two loops unless we use the improved 
lattice Laplacian, as O(a^) errors would already appear at two loop level. The two loop 
results require several more graphs and the inclusions in one loop graphs of the one loop 
mass and coupling counterterms, see Fig. Three more integrals are needed, and 
their evaluation is detailed in the appendix. The complete two loop renormahzation is 
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Figure 17: All required two loop graphs and one loop graphs with one loop counterterm 
insertions, shown as heavy dots on lines (mass coiintcrterms) or at vertices (coupling 
counterterms) . The last eight graphs cancel in pairs. Diagrams (a), (d), and (e) are 
not separately IR convergent; diagram (d) must be distributed between the other two 
to produce IR convergent integrals. 

given by 
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S{<p')-S{<p')n = (^1 + ^2)^. (1.19) 

Three new numerical constants appear here; their numerical values are Ci = .0550612, 
C2 = .0334416, and C3 = —.86147916. Sm^ and can be defined when using 

Eq. (|1.1|) , in which case C^^u = .08848010; but the other improvements do not make 
sense and should not be applied if that Laplacian is used. We can again check that the 
form of SXi and 6X2 is consistent at the two special values Ai = A2 and Ai = 3A2. 

It is possible in principle to extend the improvement scheme to O(a^), by making 
a three loop calculation. However, at this order it is necessary to include mixing of 
different dimensions of operator insertions and to include counterterms for radiatively 
induced high dimension operators in the Lagrangian. The calculation of the graphs 
also becomes significantly more challenging. The improvement presented here is suffi- 
cient for our purposes; the leading lattice spacing errors in measurables related to the 
strength of the phase transition now first appear at O(a^). 



A. 2 Two loop graphs 

Here we detail the calculation of Ci, C2, and C3. 

We begin with the two loop vertex correction. Figure ^ graph (a). The required 
integral, including the appropriate amount of the one loop counterterm graph (d), is 



167r^ Jk,BZ (P)2 yji Bz i2(j. _|_ 47r j 7fc,SR3 k'^ Ji,m P(k + /) 




(1.20) 



Here we use the shorthand in which the integration limit also lists which variable is 
being integrated over; so Jk^BZ means /[_^ j^]3((i^fc/(27r)^). The first step is to evaluate 
the inner continuum integral, which can be done by standard Feynman parameter 
methods; 

Ak3 I2{k + If ~ Jo "2^io (/2 + a(l-a)F)2 ~ 8k ' ^^''^^^ 
Then we re-arrange the original integral into three parts, 

1 1 /■ 1 / 1 1 \ /■ 1 



lk,BZ {k'^Y yji,BZ l2(^ir^l-^^ 8k An j Jk,Bz8k\(k'^f k^ J Jk,^^-Bz8k^' 

(1.22) 

The first integral is IR well behaved because the two counterterms cancel the / in- 
tegral up to a k'^ correction, which in the small k limit is .0125438fc^/47r. The inte- 
grals can be performed by quadratures using adaptive mesh refinement techniques and 
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Richardson extrapolation. The first integral gives .03 6 03/1677^ and the second gives 
.054568958/1677^. The last integral, over 3?^ — BZ^ can be re-arranged into 



1677^ JO 



dx / dy- 



1 



-.035507296027. 



1 + a;2 + j/2)5/2 



I6772 



;i.23) 



These sum to give Ci = .0550612. 

Besides this graph there is graph (e), which gives 



1 



k,BZ (k 



2\2 



1 



;i.24) 



which is not IR convergent; however, including -2 times the counterterm diagram (d), 



k,BZ {k 



2\2 



k,BZ (k 



2\2 



fe,sR3 k^ 



;i.25) 



gives — (^/477)2; no new integrals are required. It is a nontrivial check on the calculation 
that the sum of the coefficients arising from diagrams (a) and (e) precisely absorb 
diagram (d). 

The next integral is the 0{p'^) contribution from the setting sun diagram (b). 



I6772 



lim — < 



k,BZ 



(k + p) 

V 1 



1 



1 



{k + pf ky \Ji,x-i P{k + iy 



;i.26) 



The first trick is to note that 



k,BZ 



Xk + pY 



1 



:i.27) 



just by shifting the integration variable for the first term; so we may add — ^/477 to 

the term in the second parenthesis of the first line of Eq. (|1.26|) . This prevents IR 

— 2 

divergences in what follows, so we are free to expand l/{k + p) to second order in p; 
after averaging over directions for p, we find 



^ 1 ^2 



1 / 8 sin kj —sin 2fc;; \ ^ I ^ f 
3 l^i \ 3 ) 3 ^* V 



4 cos fcj— cos 2ki 
3 



2\3 



{k 



2\2 



p^M{k). (1.28) 
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The equivalent expression in the continuum case is (p^/3)//c^. Re-arranging the terms 
a httle, we can write 



Co 



167r2 



+ 



24 A,sR3_B2 

[ M{k) 
Jk — 



lk,BZ 



+ 



BZ 



k,BZ 8k 
1 

p{k'^iy 



M{k) - 

1 



(l/3) \ 
47r 



;i.29) 



The first integral is (1/3) of Eq. (|1.23|) . The second gives . 031075 7695/167?^ and the 
last gives .0142016/167r2, so C2 = .0334416. 

Next we must compute the 0{p^) part of the setting sun diagram. The continuum 
diagram is log IR and UV divergent, while the lattice diagram is only log IR divergent. 
It is convenient to IR regulate both by introducing a mass on one line. In this case 
the continuum integral can be performed in MS, leaving a lattice integral minus an 
analytically determined counterterm ^7\, ^ . Choosing to separate the renormalization 
dependence along with the same finite constant as in the previous literature [^, ^ 
the constant C3 is given by 



167r2 



lim 

m— >0 



k,BZ kP' + rn? Ji^BZ i2(j^ j^l^^ 167?^ 



1 , 6 
- + In — 

2 m 



;i.3o) 



The annoying feature of this expression is the logarithm. To remove it, we add and 
subtract l/8fc to the integral over /. The integral 



k,BZ k^ + 



LBZ 



P{k + l) 



1 

8fc 



;i.31) 



is IR convergent and the m — limit may be taken immediately. It evaluates to 
— .06858432/1671^, unless we use the unimproved lattice Laplacian, in which case it is 
.60953343/167r2. We re-arrange the remaining terms to be 



1 



k,BZ yk"^ + im? 




1 



1 



k,BZ 8k{k'^ + m?) IQ-k"^ 



6 

m. 



:i.32) 



Again, for the first integral the m ^ limit may be taken immediately, and the 
numerical value is .161799607/167?^, or .43364112015/167?^ if we use the unimproved 
Laplacian. For the last integral, we cut the integration region into the ball of radius % 
and the region within the Brillouin zone but outside the ball: 



lk,BZ 8k{k'^ 



k'^dk 



+ 



27?^ Jo 8k{k'^ + TV?) Jk,BZ 8k{k'^ + 



-e 



7r) . (1.33) 
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The former may be promptly integrated to give ln(7r/m)/167r^ plus terms power sup- 
pressed in m; when added to (— l/167r^)(ln(6/m) + 1/2) this cancels the ln(m), leaving 
(l/167r^)(ln(7r/6) — 1/2). The final integral has had the small k part of the integration 
range removed, so again the m — limit may be taken. It can then be reduced to 

1 r 1 12 /"""/^ /•arctan(sec (/)) 

-— / — ln(i?(cube) -/?(ball)) = -— — / rf0 / sin(^)c/^ln(sec(^)) 
lovr^ J 4:TT lovr^ n Jo Jo 

(1.34) 

which numerically equals .19233513195/167r^. Note that at no point have we had to 
deal numerically with an integral which is log divergent in m, or which still contains 
m at all. 

Combining terms gives C3 = —.86147916, unless we use the unimproved lattice 
Laplacian, in which case it is C3 = .08848010. In the notation of [0, C3 is called 
(. Note that, unlike Ci and C2, C3 is relatively large. Similarly, ^ is small but S is 
large. This means that the radiative 0(a) and O(a^) corrections to quantities which 
do not renormalize in the continuum are all small, but the corrections to the mass 
are larger. The size of C3 also depends on a somewhat arbitrary choice to make it 
accompany ln(6/a/i). 
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